home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1995 April / macformat-023.iso / Shareware City / Developers / RLaB / toolbox / fmin.r < prev    next >
Encoding:
Text File  |  1994-09-20  |  3.0 KB  |  124 lines  |  [TEXT/RLAB]
  1. //---------------------------------------------------------------------------
  2. //  fmin.r:
  3. //  x = fmin(f,ax,bx,tol,print)
  4. //
  5. //  find the minimum of the function f(x) in the interval [ax,bx]
  6. //
  7. //  f = The real-values function of a single variable.
  8. //
  9. //  ax,bx = left and right endpoints of search interval
  10. //
  11. //  tol = maximum length of interval of uncertainty of minimum
  12. //        tol >= 0
  13. //
  14. //  print = an optional last agrument, if nonzero a printed tracing of
  15. //          the steps is given
  16. //
  17. //  x = the value of x approximating where f attains a minimum in [ax,bx]
  18.  
  19. //  This algorithm is from: "Computer Methods for Mathematical
  20. //  Computations", Forsythe, Malcolm, and Moler, Prentice-Hall, 1976.
  21.  
  22. //  Adapted by: Duane Hanselman, University of Maine, (207)-581-2246
  23. // initialization
  24. //---------------------------------------------------------------------------
  25.  
  26. fmin = function ( f, ax, bx, tol, _print)
  27. {
  28.   global (eps)
  29.  
  30.   if (!exist (_print)) { print = 0; else print = _print; }
  31.   num = 1;
  32.   seps = sqrt(eps);
  33.   c = 0.5*(3.0 - sqrt(5.0));
  34.   a = ax; b = bx;
  35.   v = a + c*(b-a);
  36.   w = v; x = v;
  37.   d = 0.0; e = 0.0;
  38.   fx = f (x);
  39.   if (print) { fmin_data = [1, x, fx] ? }
  40.   fv = fx; fw = fx;
  41.   xm = 0.5*(a+b);
  42.   tol1 = seps*abs(x) + tol/3.0;
  43.   tol2 = 2.0*tol1;
  44.  
  45.   // main loop
  46.  
  47.   while (abs (x-xm) > (tol2 - 0.5*(b-a)))
  48.   {
  49.     num = num+1;
  50.     gs = 1;
  51.  
  52.     // is parabolic fit possible
  53.  
  54.     if (abs(e) > tol1)        // yes, so fit parabola
  55.     {
  56.       gs = 0;
  57.       r = (x-w)*(fx-fv);
  58.       q = (x-v)*(fx-fw);
  59.       p = (x-v)*q-(x-w)*r;
  60.       q = 2.0*(q-r);
  61.       if (q > 0.0) { p = -p; }
  62.       q = abs (q);
  63.       r = e;
  64.       e = d;
  65.  
  66.       // is parabola acceptable
  67.  
  68.       if ((abs (p) < abs (0.5*q*r)) && (p > q*(a-x)) && (p<q*(b-x)))
  69.       {
  70.         // yes, parabolic interpolation step
  71.         d = p/q;
  72.         u = x+d;
  73.         step = "   num        x        fx       parabolic";
  74.  
  75.         // f must not be evaluated too close to ax or bx
  76.  
  77.         if (((u-a) < tol2) || ((b-u) < tol2))
  78.         {
  79.           si = sign(xm-x) + ((xm-x) == 0);
  80.           d = tol1*si;
  81.         }
  82.       else    // not acceptable, must do a golden section step
  83.         gs=1;
  84.       }
  85.     }
  86.     if (gs)    // a golden-section step is required
  87.     {
  88.       if (x >= xm) { e = a-x; else e = b-x; }
  89.       d = c*e;
  90.       step = "   num        x        fx       golden   ";
  91.     }
  92.  
  93.     // f must not be evaluated too close to x
  94.  
  95.     si = sign(d) + (d == 0);
  96.     u = x + si * max ([abs (d), tol1]);
  97.     fu = f (u);
  98.     if (print) { fmin_data = [num, u, fu] ? disp(step) ? }
  99.  
  100.     // update a, b, v, w, x, xm, tol1, tol2
  101.     if (fu <= fx)
  102.     {
  103.       if (u >= x) { a = x; else b = x; }
  104.       v = w; fv = fw;
  105.       w = x; fw = fx;
  106.       x = u; fx = fu;
  107.     else    // fu > fx
  108.       if (u < x) { a = u; else b = u; }
  109.       if ( (fu <= fw) || (w == x) )
  110.       {
  111.         v = w; fv = fw;
  112.         w = u; fw = fu;
  113.       else if ( (fu <= fv) || (v == x) || (v == w) ) {
  114.         v = u; fv = fu;
  115.       }}
  116.     }
  117.     xm = 0.5*(a+b);
  118.     tol1 = seps*abs(x) + tol/3.0;
  119.     tol2 = 2.0*tol1;
  120.   }
  121.  
  122.   return x
  123. };
  124.